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^ : Abstract 

D , Granular elasticity, an elasticity theory useful for calculating static stress distribution in gran- 
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ular media, is generalized to the dynamic case by including the plastic contribution of the strain. 
A complete hydrodynamic theory is derived based on the hypothesis that granular medium turns 

f-j— > ■ transiently elastic when deformed. This theory includes both the true and the granular tempera- 
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c/3 , tures, and employs a free energy expression that encapsulates a full jamming phase diagram, in the 
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^ ' space spanned by pressure, shear stress, density and granular temperature. For the special case 
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I ■ of stationary granular temperatures, the derived hydrodynamic theory reduces to hypoplasticity, a 
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state-of-the-art engineering model. 
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I. INTRODUCTION 



Widespread interests in granular media were aroused among physicists a decade ago, 
stimulated in large part by review articles revealing the intriguingand improbable fact that 



13, SB- 



The resultant 



something as familiar as sand is still rather poorly understood 

collective efforts have since greatly enhanced our understanding of granular media, though 

the majority of theoretic considerations have focused either on the limit of highly excited 
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gaseous state 5|, Ig, [Zl, la, |9|], or that of the fluid-like flow [ll|, 

noteworthy and insightful simulations [iJ, ll5|, [l6| , the quasi-static, elasto- plastic motion of 



13|. Except in some 



dense granular media - of technical relevance and hence a reign of engineers - received less 
attention among physicists. 

This choice is due at least in part to the highly confusing state of engineering theo- 
ries, where innumerable continuum mechanical models compete, employing vastly different 
expressions. Although the better ones achieve considerable realism when confined to the ef- 
fects they were constructed for, these differential equations are more a rendition of complex 
empirical data, less a reflection of the underlying physics. In a forthcoming book on soil 
mechanics by Gudehus, phrases such as morass of equations and jungle of data are deemed 
apt metaphors. 



Most eng ineering theories are elasto-plastic 17 
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ones 
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19|, though there are also hypoplastic 
2l| . which manage to retain the realism while being simpler and more explicit. Both 



adhere to the continuum mechanical formalism laid down by Truesdell 22|, |23|] who, starting 
from momentum conservation, focuses on the total stress 0",^, and considers its dependence 
on the variables: strain Eij, velocity gradient V^fj and mass density p. Frequently, an 
explicit expression for aij appears impossible, incremental relations are then constructed, 
expressing dt<Jij in terms of aij, ViVj, p. Because the macroscopic energy (such as its kinetic 
or elastic contribution) dissipates, Truesdell does not include energy conservation in his 
standard prescription. 

In contrast, conservation of total energy is an essential part of the hydrodynamic approach 
to macroscopic field theories, pioneered in the context of superfiuid helium by Landau 24 1 



and Khalatnikov 



25| . The total energy w they consider depends, in addition to the relevant 



macroscopic variables such as p and fj, also on the entropy density s. (There are different 
though equivalent ways to understand s. The appropriate one here is to take it as the 



summary variable for all implicit, microscopic degrees of freedom. So the energy change 
associated with s, always written as {dw/ds)ds = Tds, is the increase of energy contained in 
these degrees of freedom - what we usually refer to as heat increase.) When the macroscopic 
energy dissipates into the microscopic degrees of freedom, the change in entropy is such that 
the increase in heat is equal to the loss of macroscopic energy, with the total energy w being 
conserved. 

y employed to account for 



The hydrodynamic approach 26|, |27| has since been successful 
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superconductors 
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3\, macroscopic electro-magnetism 4J, 
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56| . Transientlv elastic media such 
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as polymers are under active consideration at present [57, 

The main advantage of the hydrodynamic approach is its stringency. In the Truesdell 
approach, apart from objectivity, few general constraints exist for the functional dependence 
of aij or dtCTij. Therefore, one needs to rely entirely on experimental data input. In contrast, 
the structure of the hydrodynamic theory is essentially given once the set of variables is 
chosen. This is a result of the constraints provided by energy conservation, which enables one 
to fully determine the form of all fluxes, including especially the stress aij. These expressions 
are given in terms of the energy's variables and conjugate variables, they are valid irrespective 
what form the energy w has. (If w is a. function of s,p,eij, the conjugate variables are the 
respective derivative: Temperature T = dw/ds, chemical potential /i = dw/dp, and elastic 
stress vTjj = —dw/duij.) We refer to the fluxes as the structure of the theory, while taking 
the explicit form of w{s, p, Uij) as a scalar material quantity. 

There is little doubt that constructing a granular hydrodynamic theory is both useful and 
possible: Useful, because it should help to illuminate and order the complex macroscopic 
behavior of granular solid; possible, because total energy is conserved in granular media, as 
it is in any other system. When comparing agitated sand to molecular gas, it is frequently 
emphasized that the kinetic energy, although conserved in the latter system, is not in the 
former, because the grains collide inelastically. This is undoubtedly true, but it does not 
rule out the conservation of total energy, which includes especially the heat in the grains, 
and in the air (or liquid) between them. 

To actually construct the granular hydrodynamic theory, we need to start from some 
assumptions about the essence of granular physics. Our choice is specified below, and 



argued for throughout this manuscript. As we shall see, it is a guiding notion complete 
enough for the derivation of a consistent hydrodynamic theory, the presentation of which is 
the main purpose of the present manuscript. On the other hand, we are fully aware that 
only future works will show whether our assumptions are appropriate, whether the resultant 
set of partial differential equations is indeed "granular hydrodynamics." 

Granular motion may be divided into two parts, the macroscopic one arising from the 
large-scaled, smooth velocity of the medium, and the mesoscopic one from the small-scaled, 
stochastic movements of the grains. The first is as usual accounted for by the hydrodynamic 
variable of velocity, the second we shall account for by a scalar, the granular temperature Tg 
- although the analogy to molecular motion is quite imperfect: The grains do not typically 
have velocities with a Gaussian distribution, and equipartition is usually violated. All this, 
as we shall see, is quite irrelevant in the present context. 

Tg may be created by external perturbations such as tapping, or internally, by nonuni- 
form macroscopic motion such as shear - as a result of both the grains will jiggle and slide. 
Then the grains will loose contact with one another briefly, during which their individual 
deformation will partially relax. When the deformation is being diminished, so will the 
associated static stress be. This is the reason granular media can sustain static stress only 
when at rest, but looses it gradually when being tapped or sheared. And our assumption 
is, this happens similarly no matter how the grains jiggle and slide, and we may therefore 
parameterize their stochastic motion as a scalar Tg. Our guiding notion is therefore: Gran- 
ular media are transiently elastic; the elastic stress relaxes toward zero, with a rate r~^ that 
grows with Tg, most simply as r^^ ~ Tg. 

In granular statics, the grains are at rest, hence T„ = 0. With r ~ T~^ infinite, granular 



stress persists forever, displaying in essence elastic behavior |61l . |62| . I63l . |6j, |65|. When 
granular media are being sheared, because the grains move nonuniformly and Tg ^ 0, the 
stress relaxes irreversibly. This is a qualitative change from the elastic, purely reversible 
behavior of ideal solids. We believe, and have some evidence, that it is this irreversible 
relaxation that lies at the heart of plastic granular flows. If true, this insight would greatly 
simplify our understanding of granular media: Stress relaxation is an elementary process, 
while plastic flows are infamous for their complexity. 



In a recent Letter 66|, some simplified equations were derived based on the above guiding 



notion. For the special case of a stationary T„, these reproduce the basic structure of 



hypoplasticity [20], a state-of-the-art, rate-independent soil-mechanical model, and yields 
an account of granular plastic flow that is surprisingly realistic. As this agreement is a 
result of fitting merely four numbers, we may with some confidence take it as an indication 
that transient elasticity is indeed a sound starting point, from which granular hydrodynamics 
may be derived. It is not clear to us whether this starting point alone is sufficient. More 
work and exploration is needed, and especially cyclic loading, critical state, shear banding 
and tapping need to be considered. We reserve the study of these phenomena for the future. 
In this paper, we take a first step in our long march by deriving a consistent, hydrodynamic 
framework (called GSH for granular solid hydrodynamics) starting from transient elasticity. 

The paper is organized as follows. In section [III "we discuss to what extent granular media 
are elastic, or better, permanently elastic. It is well known that, although the process leading 
to a given granular state is typically predominantly plastic, the excess stress field induced 
by a small external force in a pre-stressed, static state can be described by the equations of 
elasticity. We explain why, for Tg = 0, granular elasticity in fact extends well beyond this 
limit, that it may be employed to calculate all static stresses, not only incremental ones. The 
basic reason is, without a finite Tg, there is no stress relaxation and plastic fiow. Similarly, 
if an incremental strain is small enough, producing insufficient Tg, there is too little plastic 
fiow to mar the elasticity of a stress increment. 

Then we proceed, in section IIIH to discuss jamming, a word coined to describe a system 
prevented from exploring the phase space, and confined to a single state. Although this 



idea has proven rather useful [67|, one must not forget that it is a partial view, based on 
a truncated mesoscopic model, and inappropriate for the present purpose. In this section, 
jamming is generalized and embedded in the concept of constrained equilibria. The point is, 
individual grains are unlike atoms already macroscopic. They contain innumerable internal 
degrees of freedom that are neglected in mesoscopic models |5|, la, [Zl, Isl, |9|. For instance, 
phonons contained in individual grains do explore the phase space and arrive at a distribu- 
tion appropriate for the ambient temperature. Jamming fixes only a few out of many, many 
degrees of freedom. Realizing this, the fact that grains are prevented from moving becomes 
comparable to the following textbook example: Two chambers of different pressure, sepa- 
rated by a jammed piston, and prevented from going to the lowest-energy state of equalized 
pressure. Such a system is in equilibrium and amenable to thermodynamics, albeit under 
the constraint of two constant subvolumes. Similarly, a jammed granular system at T„ = 



is also in equilibrium, not in a single state, and amenable to thermodynamics, although (as 
we shall see) under the local constraint of a given density field p{r) that cannot change even 
when nonuniform. Exploring this analogy, section IIIII arrives at a number of equilibrium 
conditions, useful both for describing granular statics and setting up granular dynamics. 

In section HVl the physics of the granular temperature is specified and developed. As 
mentioned, the energy change dw from all microscopic, implicit variables is usually subsumed 
as Tds, with s the entropy and T = dw/ds its conjugate variable. From this, we divide out 
the mesoscopic, intergranular degrees of freedom (such as the kinetic and elastic energy of 
random, small-scaled granular motion), denoting them summarily as the granular entropy 
Sg. This is necessary, because these are frequently rather more strongly agitated than the 
truly microscopic ones, Tg = dw/dsg ^ T. Note that in granular solids, we are equally 
interested in the regime Tg > T, as this is where the elasticity switches from being transient 
to permanent. In section IIVI the equilibrium condition and equation of motion for Sg are 
derived - by taking it to be an independent, macroscopic variable, without any assumptions 
about how "thermal" the associated mesoscopic degrees of freedom are. (As mentioned, 
usually they are not Gaussian and do not satisfy the equipartition theorem.) However, we 
do assume a two-step irreversibility, that the energy only goes from the macroscopic degrees 
of freedom to the mesoscopic, intergranular ones summarized in Sg, and from there to the 
microscopic, innergranular ones s. The final subsection deals with a misconception that, 
because the fiuctuation-dissipation theorem (fdt) in terms of the granular temperature 
does not usually hold, neither does the Onsager relation. The point is, the validity of fdt in 
terms of the true temperature is never in question. And the Onsager relation only depends 
on the latter. 

In section |Vl the equation of motion for the elastic strain is elucidated, and shown 
to fully determine the evolution of the plastic strain as well. In section IVIt an explicit 
expression for the free energy / is presented. This is necessary, because the energy w, or 
equivalently the free energy /, are (as discussed above) material quantities. As such, the 
free energy must be found either by careful observation of experimental data, an exercise in 
trial and error, or more systematically, through simulation and microscopic consideration. 
We proceed along the first line, making use mainly of the jamming transition that occurs 
as a function of p, Tg, Uij, to find this expression. Section fVIII presents the formal derivation 
of the hydrodynamic theory. The resulting equations are then applied to reproduce the 



hypoplastic model in section IVIIIi Finally, section IIXI gives a brief summary. 

II. SAND A TRANSIENTLY ELASTIC MEDIUM 

Granular media possess different phases that, depending on the grain's ratio of elastic 
to kinetic energy, may loosely be referred to as gaseous, liquid and solid. Moving fast and 
being free most of the time, the grains in the gaseous phase have much kinetic, but next to 
„one elastic, en.gy Q. fl fl fl 3- .. the dense. U,uid phase, say . ehute flow, the.e is 
less kinetic energy, more durable deformation, and a rich rheology that has been scrutinized 
recently [ll|, ll2|, [l3| . In granular statics, with the grains deformed but stationary, the energy 



is all elastic. This state is legitimately referred to as solid because static shear stresses are 
sustained. If a granular solid is slowly sheared, the predominant part of the energy remains 
elastic, and we shall continue to refer to it as being solid. 

When a granular solid is being compressed and sheared, the deformation of individual 
grains leads to reversible energy storage that sustains a static, elastic stress. But they also 
jiggle and slide, heating up the system irreversibly. Therefore, the macroscopic granular 
strain field Eij = Uij + Pij has two contributions, an elastic one Uij for deforming the grains, 
and a plastic one pij for the rest. The elastic energy wi{uij) is a function of Uij, not e^, and 
the elastic contribution to the stress aij is given as Tiij{uij) = —dwi/duij. With the total 
and elastic stress being equal in statics, a^j = ttjj, stress balance VjO"jj = may be closed 
Uij) , and uniquely determined employing appropriate boundary conditions. 



with TTij = TLj^ 



Our choice 
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m 



63| for the elastic energy wi = Wi{uij) is 

w, = Va (^BA' + AuA ^ B^/A Q A^ + ^ ) , (1) 



2 



TT. 



.. - -|^ = v^(SA S.^ - 2^0 + ^^'^^^' (2) 

where A = — m«, u^ = u^jU^j, u^j = Uij — ^Uu6ij. Three classical cases: silos, sand piles 
and granular sheets under a point load were solved employing these equations, producing 



rather satisfactory agreement with experiments |64J . l65l | . The elastic coefficient B, a measure 
of overall rigidity, is a function of the density p. Assuming a uniform p (hence a spatially 
constant B), the stress at the bottom of a sand pile is (as one would expect) maximal at the 
center. But a stress dip appears if an appropriate nonuniform density is assumed. Because 
the difference in the two density fields are plausibly caused by how sand is poured to form 

8 



the piles, this presents a natural resolution for the dip's history dependence, long considered 
mystifying. 

Moreover, the energy Wi is convex only for 



■Uo 



/A<v^, or 7r,/P<^/27e, (3) 



(where P = ^nu, tt^ = 7r°-7i"fj, tt^j = rcij — ^nuSij,) implying no elastic solution is stable 
outside this region. Identifying its boundary with the friction angle of 28° gives GJ, l65| 



e ^ 5/3 (4) 

for sand. Because the plastic strain pij is clearly irrelevant for the static stress, one may 
justifiably consider granular media at rest, say a sand pile, as elastic. 

If this sand pile is perturbed by periodic tapping at its base, circumstances change qual- 
itatively: Shear stresses are no longer maintained, and the conic form degrades until the 
surface becomes flat. This is because part of the grains in the pile lose contact with one 
another temporarily, during which their individual deformation decreases, implying a di- 
minishing elastic strain Uij, and correspondingly, smaller elastic energy wi{uij) and stress 
T^ij{uij)- The system is now elastic only for a transient period of time. The typical example 
for transient elasticity is of course polymer, and the reason for its elasticity being transient 
is the appreciable time it takes to disentangle polymer strands. Although the microscopic 
mechanisms are different, tapped granular media display similar macroscopic behavior, and 
share the same hydrodynamic structure. 

When being slowly sheared, or otherwise deformed, granular media behaves similarly to 
being tapped, and turn transiently elastic. This is because in addition to moving with the 
large-scale shear velocity fj, the grains also slip and jiggle, in deviation of it. Again, this 
allows temporary, partial unjamming, and leads to a relaxing Uij. 

One does not have to assume that this deviatory motion is completely random, satisfying 
equipartition and resembling molecular motion in a gas. It suffices that the elasticity turns 
transient the same way, no matter what kind of deviatory motion is present. In either cases, 
it is sensible to quantify this motion with a scalar. Referring to it as the granular entropy 
or temperature is suggestive and helpful. The granular entropy Sg thus introduced is an 
independent variable of GSH, with an equation of motion that accounts for the generation 
of Tg by shear flows, and how the energy contained in Tg leaks into heat. Only when Tg is 
large enough, of course, is granular elasticity noticeably transient. 



III. JAMMING AND GRANULAR EQUILIBRIA 

Liquid and solid equilibria are first described, then shown to correspond to the unjammed 
and jammed equilibria of granular media. 

A. Liquid Equilibrium 

In liquid, the conserved energy density w{s,p,gi) depends on the densities of entropy s, 
mass p, and momentum Qi = pvi. The dependence on Qi is universal, given simply by 

w{s,p,gi) = Wo{s,p)+g'f/{2p), (5) 

leaving the rest-frame energy wq to contain the material dependent part. Its infinitesimal 
change, dwo = {dwo/ds)ds + {dwo/dp)dp, is conventionally written as 

dwo = Tds + pdp, (6) 

by defining 

T = dwo/ds\p, p = dwo/dp\s. (7) 

It is useful to note that given Eq ([5]), the relation dw/dpls^g^ = p — v'^/2 holds, hence 

dw = Tds + {p- v^/2)dp + Vidgi. (8) 

Consider a closed system, of given volume V = Jd^r, energy j wd^r, and mass j pd^r. 
Whatever the initial conditions, it will eventually arrive at equilibrium, in which the entropy 
J sd'^r is maximal, or equivalently, at minimal energy for given entropy, mass and volume. 
To obtain the mathematical expression for this final state, one varies J wd^r for given J sd^r 
and j pd^r, arriving at the following equilibrium conditions, 

V,T = 0, Vi/i = 0. (9) 

Being expressions for optimal distribution of entropy and mass, these two conditions may 
respectively be referred to as the thermal and chemical one. 

In mathematics, Eqs iQ are referred to as the Euler-Lagrange equations of the calculus 



of variation. The calculation is given in Appendix |Al More details may be found in 68| , 



10 



in which three additional conserved quantities: momentum J gid'^r, angular momentum 
j {r X g)id^r, and booster /(prj — git)d^r were also considered, adding a motional condition, 



V. 



«j 



(ViVj + VjVi)/2 = 0, (10) 



and altering the chemical one to dtVi + Vi(/i — f^/2) = 0. We focus on Eqs iQ here. 

Including gravitation, the energy is wq = wq + 0, with Gk = — Vi0 the gravitational 
constant pointing downwards. The generalized chemical potential is 

/i(p) = dwo/dp = /i + (/), (11) 

while chemical equilibrium, Vj/i = 0, is 

Vi/i = Gi. (12) 

This implies a nonuniform density represents the optimal mass distribution minimizing the 
energy (or maximizing the entropy). With the pressure given as Pt = —Wq + TS + fip, see 
Appendix |Al the condition for mechanical equilibrium, 

ViPr = sViT + pVifi = pGi (13) 

is a combination of the thermal and chemical ones. 

B. Solid Equilibrium 

In solids, if the subtle effect of mass defects is neglected, density is not an independent 
variable and varies with the strain (for small strains) as 

dp/p = -duu. (14) 

Defining ttj^ = —dwo/duij\s, we write the change of the energy as 

dwQ{s,Uij) = Tds — TTijduij. (15) 

Maximal entropy, with the displacement vanishing at the system's surface, implies the fol- 
lowing thermal and mechanical equilibrium conditions (see Appendix [Al). 

ViT = 0, V.TT,, = 0. (16) 

11 



So force balance is, in the complete world including the innergranular degrees of freedom, an 
expression of maximal entropy - quite analogous to uniformity of temperature. It implies 
the overwhelming dominance of phonon distribution that satisfies force balance, and the 
rarity of phonon fluctuations that violate it. 

Including gravitation, the total energy is given as dwo{s,Uij) = Tds — nijduij, with 
TTjj = TTij + p(f), and mechanical equilibrium becomes 

VjTTij = pGi (17) 

C. Granular Equilibria 

Depending on whether Tg is zero or finite, sand flip-flops between the above two types 
of behavior. The density is an independent variable, because the grains may be differently 
packaged, leading to a density variation of between 10 and 20% at vanishing deformation. 
So the energy depends on all three variables, 

dwo(s, p, Uij) = Tds + pdp — iiijduij. (18) 

If Tg is finite, the elastic stress ttjj relaxes until it vanishes. The equilibrium conditions are 
therefore, including gravitation, 

V.T = 0, ViPT = pG,, vTi, = 0, (19) 

similar to that of a liquid, with ViPr = pGi (or Vj/i = Gi) enforcing an appropriate density 
field, and tt^ = forbidding any free surface other than horizontal. 

For vanishing Tg, sand is jammed, implying two points: First, ttjj no longer relaxes; 
second, without slipping and jiggling, the packaging density cannot change, and the density 
is again a dependent variable, dp/p = — dn«. The suitable equilibrium conditions, as derived 
in Appendix El are 

ViT = 0, VjiPrSij + TCij) = pG^, (20) 

which allow static shear stresses and tilted free surfaces. So, although jammed states are 
prevented from arriving at the liquid-like conditions of Eqs flT9|) . they do possess reachable 
thermal and mechanical equilibria. 
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If the energy (as given in section lIVp depends in addition on the granular entropy, dw 
Tds + Tgdsg + ■ ■ ■ , the pressure contribution Pt (see section IVII A|) is 



Pt = -Wo + Ts + TgSg + fip = -f + fip, (21) 

with ViPr = sViT + SgViTg + pVifi. (22) 

IV. GRANULAR TEMPERATURE Tg 

Granular temperature is not a new concept. Haff, at the same time Jenkins and Savage p, 
la, LZI, la, |9|] , introduced it in the context of granular gas, taking (in an analogy to ideal gas) 
Tg ~ Wkin, where U7kin is the kinetic energy density of the grains in a quiescent granular 
gas. With Tg = dwi^i^/dsg ~ dTg/dsg, the granular entropy is Sg ~ InT^,. As discussed 
above, granular temperature is also a crucial variable in granular solids. But one must not 
expect this gas-like behavior to extend to the vicinity of Tg = T: As the system, if left 
alone, always returns to Tg = T, the energy must have a minimum there. And something 
like ti; ~ s^ ~ (Tg — T)^ and Sg ^ Tg — T would be more appropriate. (Neither for ideal 
gases does Sg ~ InT^ persist for all temperature. Excluding a phase transition, quantum 
effects become important before T = is reached.) 

A. The Equilibrium Condition for Tg 

The energy change dw from all microscopic, implicit variables is generally subsumed as 
Tds, with s the entropy and T = dw^/ds its conjugate variable. From this, we divide out 
the intergranular energy of the random motion of the grains, denoting it as TgdSg, 

dWO = Td(S - Sg) + TgdSg = T dS + (Tg - T)dSg. (23) 

The first expression distinguishes between two heat pools: s — Sg and Sg, with the latter 
rather more strongly excited, Tg 3> T. The second expression, algebraically identical, takes 
w as a function of s and Sg, with Tds being the total heat if all degrees were at T, and 
(Tg — T)dsg the increase in energy when some of the degrees are at Tg. If unperturbed, a 
stable system will always return to equilibrium, at which the second pool is empty, Sg = 0. 
This implies the free energy f = wq — Ts has a minimum at Sg = 0. Assuming analyticity. 
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we expand the free energy f(T, Sg) around Sg = 0, arriving at 

f = fo{T) + sl/{2bp), (24) 

where 6 is a positive material parameter, a function of p and Uij. [The factor p will turn out 
later to be convenient.] With d/ = —sdT + (Tg — T)dsg we have 

Tg=Tg-T = df/dSg\T = S g / {h p) (25) 

that vanishes in equilibrium 

fg^Tg-T = 0. (26) 

We shall employ the Legendre transformed potential, f(T, Tg) = f(T, Sg)—TgSg, below (that 
has a maximum rather than a minimum at Tg = T), 

fiT,fg) = foiT)-bpf;/2. (27) 

Because an improbably high Tg is implied by any random motion of the grains, neglecting 
T in comparison to Tg or taking Tg ^ Tg is frequently a good approximation, though not 
close to Tg = 0. So it is prudent not to implement it while deriving the equations. 

B. The Equation of Motion for Sg 

Being a macroscopic, non-hydrodynamic variable, Sg must first of all obey a relaxation 
equation, —dtSg = -ydf /dsg = -yTg. Since this relaxation is typically slow, Sg also displays 
characteristics of a quasi-conserved quantity, and removal of local accumulations is accounted 
for by a convective and a diffusive term, 

-dtSg = Vi[SgVi - KgVifg] + 'jTg (28) 



V^{sgVi) + (1 - xV')sg/r, 



9' 



where Tg = bp/'j is the relaxation time, while x — ^/^g/l is the characteristic length 
associated with the diffusion. (The second line of Eq fl25]) assumes Kg, 7 = constant.) If Tg 
is held at Tq at the boundary x = 0, and allowed to relax for a; > 0, the field Sg ~ Tg{x) 
obeys (1 — x^V^)sc, = in the stationary limit dtSg,Vi = 0, and decays as 

Tgix)=T + Toexp{-x/x). (29) 

14 



Eq (l28l) is not complete. To see this, consider first the true entropy s. In hquid, s is 
governed by a balance equation with a positive source term R that is fed by shear and 



compressional flows, and by temperature gradients [24 1 



dfS + V^isV^ - KWiT) = R/T, (30) 

^ = ^<4 + c4 + «(v.r)2, (31) 

where w°- is the traceless part of Vij = ^{ViVj + V^f ») and vu its trace; 77, ^ > are the shear 
and compressional viscosity, respectively, and k > the heat diffusion coefficient. Entropy 
production R must vanish in equilibrium and be positive definite off it. The thermodynamic 
forces VjT and Vij also vanish in equilibrium [see Eqs (PITOJ) ]: off it, they may be taken 
to quantify the "distance from equilibrium." The entropy production R increases with this 
distance and may be expanded in VjT and Vij. The given terms are the lowest order, positive 
ones that are compatible with isotropy. 

In granular media, equilibrium conditions are more numerous than in liquid. As discussed 
in section [Til Ct these are, in addition, the vanishing of ttj^, VjTTjj, and Tg, hence we have 

R = Vv%vl + C4 + <'^^Tf + if'g (32) 

+/3«)^ + A4 + /3^(V,7r,,)^. 

Three additional points: (1) Being an expansion in the thermodynamic forces, the transport 
coefficients ri,(, t^, t^gyl, P, Pi, P^ may still depend on the variables of the energy, T,Tg,p, 
TTfj and TT^ = ^Ij'^ij' ^^^ ^'-'^ "-"^ ^^^ forces themselves, such as VjT or Vij. (2) More terms 
are conceivable in Eq (!32l) . say aiV iTV jiiij or kitTjjVjTVjT. These may be included when 
necessary. (3) The above reasoning leaves the question open why Vj/i does not contribute 
to R, not even in liquid - or more precisely, why the coefficient preceding (Vj/i)^ always 



vanishes. The answer is given in 68|, though there have been some recent controversies 
about it, see |69| and references therein. 

The granular entropy Sg should obey a balance equation with the same structure, 

dtSg + Vi{SgVi - KgWifg) = Rg / Tg , (33) 

though the source term Rg has positive as well as negative contributions: Two positive ones 
from shear and compressional flows, and the negative relaxation term discussed in Eq ( l28ll . 

R9 = VgVlvl + CgVl + Kg{V,fgf - ^f^. (34) 
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The fact that the coefficient preceding Tg is 7 both in Eq (l32l) and (1M|) derives from energy 
conservation: Taking the system to be uniform, we have dtw = TdtS+TgdtSg = R+Tg{—jTg) . 
So dtW = imphes R = 'jTg. It expresses the fact that the same amount of heat leaving Sg 
must arrive at s. A direct consequence for the stationary case, Rg = 0, is 

^f^ = rigv^^v^^ + C,4, (35) 

quantifying how much Tg = Tg — T is excited by shear or compressional flows. 

In dry sand, the granular viscosities rjg, Q probably dominate, while 7], ( are insignificant 
- though the latter should be quite a bit larger in sand saturated with water: A macroscopic 
shear flow of water implies much stronger microscopic ones in the fluid layers between the 
grains, and the energy dissipated there goes to the true entropy s, instead of to Sg ffist. 

C. Two Fluctuation-Dissipation Theorems 

There are many in the granular community who dispute the validity of the Onsager 
reciprocity relation in granular media, enlisting any of the following three reasons: (1) The 
fluctuation-dissipation theorem (fdt) does not hold. (2) The microscopic dynamics is not 
reversible. (3) Sand is too far off equilibrium. 

Careful scrutiny shows that none of these arguments holds water. First, with F denoting 
the free energy, fluctuations say of the volume are always given as 

{AV^) = T{d^F/dV^)-^ = T{-dP/dV)-\ (36) 

Jammed sand, similar to a copper block, undergoes volume fluctuations as described by 
Eq ( l36l) . When sand is unjammed, Eq ( l36l) still holds, though F now depends on T^, such 
as given in section IVll In granular media, T is frequently replaced by T^, 

{AV^) = Tg{-dP/dV)-\ (37) 

This "fdt" is indeed highly questionable, because Tg frequently behaves rather differently 
from the true temperature. However, the crucial point here is, the validity of the Onsager 
relation depends on Eq (1361) . not Eq ( 1371) . 

Second, the dynamics typically employed in granular simulations is indeed irreversible, 
but only as a result of a model-dependent approximation that treats grains as elementary 
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constituent entities. The true microscopic dynamics that resolves the atomic building blocks 
of the grains remains reversible. And this is the basis for the Onsager relation. 

Third, "too far off equilibrium" is not convincing, as turbulent fluids, truly far off equi- 
librium, are known to obey them. Some argue that sand, whether jammed or in motion, 
are always far from equilibrium. Yet as the careful discussion in section IIIII shows, this is 
an inappropriate view. Granular media are not always far from equilibrium, they just have 
different ones to go to - solid-like if jammed and liquid-like if unjammed. 

V. ELASTIC AND PLASTIC STRAIN 

As discussed in section [TTl the elastic strain Uij accounts for the deformation of individual 
grains, while their rolling and sliding is described by the plastic strain pij. Together, they 
form the total strain 

Eij = Uij+pij. (38) 

The elastic energy w{uij) is a function of Uj^, not of e^j, and the elastic stress is given as 
TTij{uij) = —dw/duij. When Tg is finite, the elastic strain relaxes, 

dtUij - Vij = -Uij/r. (39) 

implying a diminishing elastic strain Uij, and correspondingly, smaller elastic energy w{uij) 
and stress TTij{uij). Note because the total strain is a purely kinematic quantity, dtSij = Vij, 
the evolution of the plastic strain pij is also fixed, dtPij = Vij — dtUij. 

It is the relaxation term —Uij/r that gives rise to plasticity. To see how it works, take a 
constant r and consider the following scenario. If a transiently elastic medium is deformed 
quickly enough by an external force, leaving little time for relaxation, /(mjj/t) dt ^ 0, 
we have Sij = Uij and Pij = right after the deformation. If released at this point, the 
system would snap back toward its initial state, as prescribed by momentum conservation, 
dt {pvi) + ^j'n'ij = 0, displaying thus a behavior that is clearly reversible and elastic. But if 
we hold the system still for long enough, Vij = 0, hence dt Eij = 0, the elastic part Uij will 
relax, dtUij = —Uij/r, while the plastic part grows accordingly, dtPij = —dtUij. When Uij 
vanishes, the plastic part will have completely replaced it, pij = Eij. With the elastic stress 
TTjj and energy w{uij) also gone, momentum conservation reads dt (pvi) = 0. The system 
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now stays where it is when released, and no longer strive to return to its original position. 
This is obviously what we mean by a plastic deformation. 

Next take r ~ T^^^. As discussed in the introduction, this should be appropriate for 
granular media. Assuming (for simplicity) a stationary granular temperature, or T^ = 
{j]gl'^)vijVij = (^c,/7)||'ys|p, see Eq ( 15^ . we obtain from Eq (15^ the equation. 



dtUij - f ij ~ I |f s 1 1 {-Uij) \j'r]gl-i , (40) 



the rate- independent structure of which closely resembles the hypoplastic one [201]. As a 
result, both the elastic strain Uij and the stress Oij will display incremental nonlinearity, 
ie., behave differently depending whether the load is being increased {vij > 0, \\vs\\ > 0) 
or decreased {vij < 0, ||fs|| > 0). Not surprisingly, this equation leads to plastic flows 
very similar to the hypoplastic results. However, under cyclic loading of small amplitudes, 
because Tg never has time to grow to its stationary value, the plastic term Uij/r ~ TgUij 
remains small, and the system's behavior is rather more elastic. 

The equation of motion for the elastic strain [cf. the derivation leading to Eq (177|1 ] is in 
fact somewhat more complicated and given as 

dtUij - (1 - a)vij - Xij 
= -[{uikVjVk + Viyj/2) + {i^ j)], (41) 

where d^ = dt + VkVk, and {i ^^ j) signifies the same expressions as in the preceding bracket, 
only with the indices i and j exchanged. In this equation, the term (uikVjVk) + {i ^^ j), 
important for larg e strain field and frequently negligible for hard grains, is of geometric 
origin, see [57|, l58|, |59|, l60|| for explanations. The dissipative fluxes Xij = —Uij/r — avij and 
Vi ~ ^j'^ij ^^^^ tie derived in section IVII A[ The second term is quite similar to the diffusive 



heat current kVjT, which aims to reduce temperature gradients and establish VjT = 0. 
We can take yi to be a current that aims to reduce VjTfij and establish the equilibrium 
condition, Vjttjj = 0, of Eq fl20|) . Given Eq fHTj) and dtSij + [(^jfcVjffc) + {i <-^ j)] — Vij = 0, 
the evolution for the plastic strain is again fixed. 

VI. THE GRANULAR FREE ENERGY 

As explained in the Introduction, the structure of the hydrodynamic theory is determined 
by general principles, especially energy and momentum conservation, but the explicit form 
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of the energy w is not. Although w does possess features that it must always satisfy, most 
of its functional dependence reflects the specific behavior of the material. To arrive at an 
expression for the energy of granular media, there are two obvious methods, either a mi- 
croscopic derivation, possibly via simulation, or more pragmatically, examining constraints 
from key experiments, opting for simplicity whenever possible, as we do here. 

Because we are interested in the hmit of small Tg and m^-, see Eq ^ and (12^ . and because 
the dependence on the true temperature is usually irrelevant, the difficult part is the density 
dependence of the energy. Fortunately, quite a number of known features may be used as 
input. First, there are two characteristic granular densities, the minimal and maximal ones, 
Pip and Pcpi respectively referred to as random loosest and closest packing. In the first case, 
the grains necessarily loose contact with one another when the density is further decreased; 



in the second, the density can no longer 
system is orders of magnitude stiffer 



3e increased without compression, at which point the 
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7l| . Then there is the jamming transition of 



sand, especially the so-called virgin consolidation line., which we believe is the limit beyond 
which no stable elastic solutions are possible, see Fig[T]-(a). These in conjunction with the 
density dependence of sound velocity and the pressure exerted by agitated grains contain 
sufficient information to fix the expression for the energy. 

Instead of the energy, we consider the potential f(T,Tg,p,Uij) = wo — Ts — TgSg, see 
Eq fl271) . Referring to it for simplicity also as the free energy density, we write 

/ = /o(T, p) + /i(p, u,,) + Up, fg), (42) 

h ^w,= B^/K (2AV5 + ul/i), (43) 

h = pho{l - plp,p)\-fll2), < a « 1, (44) 

where /o(T, p) is the free energy at vanishing granular temperature and elastic deformation, 
Tg, Uij = 0, while wiiuij) and f2(Tg) are the respective lowest order term. (It is a simplifying 
assumption that the temperature T enters the free energy only via /o, and not Wi, f2- This 
neglects effects such as thermal expansion that, however, may be added when necessary.) 

Being cohesionless, the grains possess no interaction energy, /o(T, p) is therefore the sum 
of the free energy in each of the grains, 

/o(r,p) = (Fi(r)/m)p, (45) 

where Fi is the free energy of a single grain, m its mass, and {Fi{T)/m) the free energy per 
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unit mass, averaged over a number of grains. 

It is important to realize that the equihbrium stress is given, once one knows what the 
free energy density f = F/V is (see Appendix R|l . 



aij = PrSij + vr, 



ij I "ij 



dif/p) 



d{l/p) 
The first term is the local expression for the more familiar one, 

d{f/p) 



'^^-^- (^^) 



9F d{fV/M) 



dV d{V/M) 



M 



d{l/p) 



= pdf/dp -f = pfi + Ts + TgSg - w. (47) 

In liquids, only this term exists, since / does not depend on Uif^ in ideal crystals, only the 
second term exists, because the density is not an independent variable, see the discussion 
in section IIIII In granular media, both terms coexist. Given the free energy / = ^ /« of 
Eq (l42l) . each term yields the pressure contribution, 

P, ^ pidfjdp) - /„ (48) 

with PT = EP^ and Po = pdfo/dp - /o = 0. 

A. The Elastic Energy 

The elastic part of the free energy, Eq (H3|) . has previously been successfully tested under 
varying circumstances, cf. the discussion in section [TTl below Eq ([2]). It is not analytic in 
the elastic strain, but does contain the lowest order terms. As it takes some deliberation to 
arrive at its density dependence and the terms of higher order in Uij, we consider them in 
two separate sections below. 

First, a conceptual point. We take any yield surface as the divide between two regions: 
One in which stable elastic solutions are possible, the other in which they are not - so a 
system under stress must flow and cannot come to rest here. Accepting this, the natural 
approach is to have a convex elastic energy turn concave at the yield surface. The idea 
behind it is, the energy is an extremum if the equilibrium conditions of section UTTl including 
especially Eq (!20|) . are met. Convexity implies the energy is at a minimum there, and 
concavity that it is at a maximum. Where Wi is concave, any elastic solution satisfying 
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Eq (120!) has maximal energy, and is eager to get rid of it. It is not stable because infinitesimal 
perturbations suffice to destroy it. 

As discussed in section HIl for B, ^= constant, Wi is convex for tt^/P < a/2/^ and concave 
otherwise, and already possesses the right form to account for the Coulomb yield line, see 
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FIG. 1: Granular yield surface, or jamming phase diagram, for Tg = 0, as a function of the 
pressure P, shear stress 0"^, and void ratio e = pa/p — 1- All thick solid lines are calculated using 
Eqs (J43l51l54p . (a): Maximal void ratio e versus pressure P, or the virgin consolidation line. The 
dotted line is an empirical formula, e = 0.679 — 0.0971n(P/0.5), with P in Mpa. The thin line 
(designated as simple model) renders Eq (j53|) . The circle at the top is the random loosest packing 
value for e. (b): The straight Coulomb yield line bends over depending on e, a behavior usually 
accounted for by the cap model in elasto-plastic theories, (c): The 3D combination of (a) and (b). 
Values for the calculation are : Bq = 7000 Mpa, pL = 0.445/9G, Pep = 0.645pG) ^i = 10~^, and 
ki = 10-5 mVkg, k2 = 1000, k^ = 0.01. 
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Fig[l]-(b). Our task now is to appropriately generalize it such that the density p is included 
as a third variable. Instead of p, the void ratio, e = pc/ P — 1, is frequently employed. It 
remains constant at elastic compressions and accounts for granular packaging only, [pc the 
bulk density of granular material, typically around 2700 kg/m'^ for sand.) 

B. Density Dependence of B 

We shall take B as density dependent, but not ^: Since the Coulomb yield line is approx- 
imately independent of the density, so must the coefficient ^ be, see Eq (jlj) . Granular sound 



velocity was measured by Hardin and Richart 72|, who found it linear in the void ratio, 



c ~ 2.17 — e. Given Eq (H3l) . the velocity of sound is c ~ yS/p, implying 

B = e„(3.17 - Pg/p)Hp/pg)- (49) 

Since this expression properly accourrts for the measured fl density dependence of the 
compliance tensor M^^i, the dependence of ;B on p seems settled |7J]. It is not, because the 
resultant wi is concave in the variables p and A, and could not possibly sustain any static 
solution. Inserting Eq fllU]) into fH31) . we find the energy violating the stability condition, 

d^B'^/yOp^ < 0, (50) 

obtained from inserting Eq (143!) with Ug = into {d'^Wi/dp^) {d'^Wi/dA'^) > (d'^wi/dpdA)'^. 
Clearly, the widely employed Hardin- Richart relation, c ~ 2.17 — e, is not accurate enough 
for a direct input into the energy. It works fine as long as the sand is jammed, T„ = 0, 



and p is only a given parameter, not a free variable - such as in the experiments of [73[], or 
when determining static stress distributions. But if a finite Tg frees the density to become 
a variable, this instability will wreck havoc with the hydrodynamic theory. We need to 
reconstruct the density dependence of B, such that the energy wi 

1. vanishes for densities smaller than the random loosest packing value (around the void 
ratio of e^p ~ 0.8 for sand of uniform grain size), or p < p^^; 

2. (as a simplification) diverges at p = pcp, the random closest packing value (around 
(icp ~ a. Do)] 

3. is convex and reproduces the Hardin- Richart relation between pip and pcp. 
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Alas, these points are more easily stated than combined in an energy expression, and no 
continuous B seems feasible: If analytic, B would be proportional to p — pep close to pip. 
More generally, we may take B ^ {p ~ Pep)"", with a positive. But the resulting energy, 
w ^ {p — p^p)"A^'^, remains concave. Only when including the divergence at pcp by taking 
B r^ (p — Pep)°'/{pcp — pY does the energy turn convex, between pcp and a density larger 
than Pip. We therefore propose 

B = B^ {P.Zl!2^ X C, for p > p,p- (51) 



Pcp- P 



i3 = 0, 



for p < Pip. 



(52) 



With an appropriate p^^ < pip, this expression renders the energy divergent at pcp, stable and 
convex up to pgp, and approximates the Hardin- Richart relation between them, see Fig. ([2]). 
(Take C = 1 for now, until it is specified otherwise in the next section.) 

C. Higher-Order Strain Terms 

Next, we consider the unjamming transition in connection with compaction by pressure 
increase, the fact that denser sand can sustain more compression before getting unjammed. 




FIG. 2: Equation (149 jl . obtained by employing the Hardin- Richart relation directly, violates the 
stability condition Eq (j50p . because d'^B"'^'^ /dp^ > for all density values. Although numerically 
similar, see insert, the expression from Eq (j5ip suitably becomes concave at pfc, and satisfies the 
stability condition between pic and pcp- The plots are calculated with />|^ = OAAbpn^Ppc = 0.645pG 



(implying pip = O.bbbpc), and Bq = 7000 Mpa, appropriate for Ham River sand [73]. 
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before elastic solutions become unstable: See the dotted line of Fig[T]-(a), depicting a well- 
known empirical formula from soil mechanics |17l . llSl . Il9| , e = eo — A In P. Referred to as 
the virgin (or primary) consolidation line, it represents the boundary that sand (at rest) 
will not cross when compressed. Instead, it will collapse, becoming more compact, with a 
smaller e, close to or at the curve, but not beyond. (Note the dotted line does not appear 
to cut the e-axis, as it should at pip - this is where sand becomes instable for any pressure. 
The discrepancy may derive from difficulties of making reliable measurements close to pep.) 
This behavior is a natural consequence of higher-order strain terms such as the next order 
ones (Ci,C2 > 0), 

-iCiA' + C^Aul), (53) 

which need to be added to wi as given by Eqs ( l43f5T]) . Consider first pure compression, 
Mg = 0. For small A, the term — (^lA^ is negligible, and Wi remains convex. But if A is 
large enough, its negative second derivative will turn wi concave, making any elastic solution 
impossible. The value of A at which this happens, grows with B - a larger third-order term 
is needed for a larger B. Now, B is smallest at p = pip, grows monotonically with p, and 
diverges at pcp- As a result, the instability line cuts the e-axis at pip, veers towards larger 
A (or larger P) at higher density , and heads for infinity at pcp, see the thin line depicted 
as "simple model" in Fig[T]-(a), drawn with a constant d = 24500 MPa. (It is of course 
possible, employing a density- dependent (i, to improve the agreement to the dotted line.) 
In Fig[l]-(b), the point of maximal pressure for a given void ratio e is located at where the 
P-axis is being cut by the associated curve. If the term ~ Au^ did not exist, these curves 
would be vertical lines. The presence of ~ Au^ reduces the value of A (or P) for growing 
Us (or as), bending the lines to the left. 



Although qua . 



in textbooks 
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itative figures of these curves that are frequently referred to as caps abound 
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19| , we did not find enough quantitative data, especially not a generally 
accepted empirical expression, that we could have compared our results to. Presumably, it 
is not easy to observe caps in dry sand. Given this lack of reliable data, we decided against 
the expansion, Eq ( l53l) . and opted for a fiexible "cap function," C of Eq (15T1) . capable of 
accounting for any possible cap- like unjamming transitions, 

2C = 1 + tanh[(Ao - A)/Ai], where (54) 

Ao = hp - k2u1 - h = k[/{e + 1) - k2u1 - k^. 
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With C ~ 1 for A ^ Ao, and C ~ for A :^ Aq, the cap function is constructed to be 
relevant only in a narrow neighborhood around Aq, for |A — Ao| < Ai ^ 10~"^, such that 
the energy's convexity is destroyed around Aq. Taking ki,k2, k^ as constant, Aq grows with 
the density and falls with u^, giving rise to the typical appearance reproduced in Fig[TJ 

Together, Eqs fl43f51f54l) give the energy density Wi, appropriate for cohesionless granular 
materials at Tg = 0. There are two contributions to the pressure, P = Pi + Pa, where 
Pi = p(dwi/dp) — wi from Eq (jUj), and TCij = —dwi/duij = -Pa% — o-gu'^j/ug. Because 
we still take A to be a small quantity. Pi ~ A^'^ may be neglected. (Similarly, terms such 
as TiikUjk ~ A^'^ from Eq ( l66l) below are also neghgible.) So the stress is simply iiik, with 
pressure and shear stress given as 

Pa = BVA{A + ^ul/A) - wiC*/Ai, (55) 

By/Aus - 2k2UsWiC* / Ai, (56) 



a. = I 



where C* = 1 — tanh[(Ao — A)/Ai], hence C* -^ away from the cap. (The terms of higher 
order in A are kept in C*, because Ai is small. This is how we make C a function relevant 
for A ^ Ao, not A -^ 0.) 

Stability is given only if the energy wi is convex with respect to its seven variables, 
p,A,u^j. As linear transformations do not alter the convexity property of any function, 
we may take the energy as W7(p, A,xi_5) where xi = ^/2uxy, x^ = \/2uxz, x^ = \/2u 



yzi 



X4 = {uxx — Uzz)/V^, 3^5 = (uxx — '2uyy + Uzz)/VQ- The characteristic polynomial Nj of the 
Hessian matrix of wj is A^y = (A — u'^^dwi/dusYNs, with N3 the characteristic polynomial 
oi wi{A,Us, p). Since u~^dwi/dus is always positive, it is sufficient to consider wi{A,Us, p). 
Requiring A^3 to have only positive eigenvalues defines the stable region in the strain space, 
spanned by A, Us,e. Using Eqs (I55ll56p . we may convert this into one in the stress space, 
spanned by P, cTs, e. The result, obtained numerically, is the yield surface plotted in Figdl 

D. Pressure Contribution Prom Agitated Grains 

Agitated grains are known to exert a pressure in granular liquid. Using the model of 
ideal gas (better: non-interacting atoms with excluded volumes), with W2 ~ pTg denoting 
the energy density of agitated grains, the pressure expression, 

Pt(p,T,)~^2/(1-p/Pcp), (57) 
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was employed and found to account realistically for the behavior o: 
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[granular liquid sand- 
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78| 



wiched between two cylinders rotating at different velocities |75l . 

In ideal gas, both the energy density w and pressure P are proportional to the temperature 
T. As a consequence, the entropy is s ~ InT, and diverges for T ^ 0. (The free energy has a 
contribution ~ TlnT that vanishes for T -^ 0.) As quantum effects become important long 
before T vanishes, the unphysical feature of a diverging entropy is inconsequential for ideal 
gases. Yet this would be a highly relevant defect for granular solids, for which important 
physics occurs at or around Tg = 0. This is the reason ideal gas is not an appropriate model 
for granular solids. The considerations of section HVl show that W2, J2 ~ ^g close to T^ = 
- implying a pressure contribution, P2 = p{df2/dp) — /2 ~ T^, see Eq ( H51) . Note first that 
P2 ~ W2 is retained, and second that because Pq = 0, Pi ~ 0, we have P-r = X] -P* ~ -^2- 

Unfortunately, the density dependence of Eq fl57j) also poses a problem, as it implies 



a free energy /2 = bQp\n{l — p/pcp){—T^/2) and a granular entropy, Sg = —df2/dTg = 
6opln(l — p/pcp)Tg, both diverging for p — > pcp. We therefore take /2 to be given as in 
Eq (H4l) . with a positive but small a. The resulting entropy is physically acceptable, and the 
pressure is easily rendered numerically indistinguishable from Eq ( l57|) . 

p apb(,Tg 

As the total pressure is now P = Pt + Pa, cf. Eq ( 155|1 . the jamming transition discussed 
above is modified. For instance, the yield condition of Eq ([3]), with ^ = 5/3, now reads 

T^s TTs /6 , ^ 

-^ = — < \ -, (60) 

Pa P-Pt~\5' ^ ^ 

implying a smaller maximal tTs for given P. On the other hand, the maximal value for the 
void ratio e is larger when P^ is present: Any given e has a maximal elastic compression A 
that will not sustain a larger e. But if P is fixed and Tg is finite, the elastic compression A 
will be appropriately smaller to sustain a larger e. This behavior is depicted in Fig El 

The jamming transition, from elastic solid to liquid, is of course no longer completely 
sharp at a finite Tg, because Tg turns the elastic body into a transiently elastic one for all 
values of stress and density. Nevertheless, there is a huge quantitative difference between 
catastrophic unjamming and the gradual process of stress relaxation. A sand pile may slowly 
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FIG. 3: Jamming transition as a function of e,as and Pt ~ T?, for Pa = 0.4 MPa. Values of 
model parameters are the same as those in FIG. 1. 



degrade, relaxing toward the fiat surface. But when turning on Tg violates Eq ( 1601) . sudden 
events such as liquefaction happen. {Pt may be substituted by the pore pressure to account 
for a similar collapse, if the soil is filled with water.) The frequently reported phenomenon 
of a primary earthquake emitting elastic waves that trigger earthquakes elsewhere [79(], may 
well be connected to Eq ( l60l) : Tg as given by Eq ( l35l) accompanies elastic waves. It may be 
sufficiently large to violate Eq (1601) if stability was precarious. 



E. The Edwards Entropy 

It is useful, with the free energy obtained in this chapter in mind, to revisit the starting 
points of Granular Statistical Mechanics (gsm), especially the Edwards entropy 80|. Taking 
the entropy S(W, V) as a function of the energy W and volume V, or dS = {l/T)dW + 
{P/T)dV, it argues that a mechanically stable agglomerate of infinitely rigid grains at rest 
has, irrespective of its volume, vanishing energy, W = 0, dW = 0. The physics is clear: 
However we package these rigid grains that neither attract nor repel each other, the energy 
remains zero. Therefore, dS = {P/T)dV, or dV = {T/P)dS = XdS. This is the starting 
expression of GSM, and X is considered the relevant quantity characterizing granular media 
at rest. The entropy 5* is obtained by counting the number of possibilities to package grains 
for a given volume, taking it to be e"^. And because a stable agglomerate is stuck in one 
single configuration, some tapping (or a similar disturbance) is taken to be needed to enable 
the system to explore the phase space. 
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In GSH, the grains are neither infinitely rigid, nor generally at rest. An elastic and a 
Tg-dependent energy contribution, denoted respectively as /i and /2, see Eq f l42l) . account 
for them. GSH also possesses a T^-switch that determines whether the system's behavior is 
solid- or liquid-like. This is clearly the generalization of phase space exploration enabled by 
tapping. That grains neither attract nor repel each other is accounted for by the stress van- 
ishing for Tg, Uij — i> 0: In this limit, in which /i, /2 = and only /o ~ p finite, there is no term 
in the energy that depends nonlinearly on the density p, hence aij = d{fo/p)/d{l/p)6ij = 0. 

Given this comparison, it is natural to ask whether GSM is a legitimate limit of GSH. The 
answer is probably no, as both appear conceptually at odds in two points, the first more 
direct, the second quite fundamental: (1) Because of the Hertz-like contact between grains, 
very little material is being deformed at first contact, and the compressibility diverges at 
vanishing compression. This is a geometric fact independent of how rigid the bulk material 
is. Infinite rigidity is therefore not a realistic limit for sand. (2) In considering the entropy, 
one must not forget that the number of possibilities to package grains for a given volume is 
vastly overwhelmed by the much more numerous configurations of the inner granular degrees 
of freedom. Maximal entropy S for given energy therefore realistically implies minimal 
macroscopic energy, such that a maximally possible amount of energy is in S (or heat), 
equally distributed among the numerous inner granular degrees of freedom. Maximal number 
of possibilities to package grains for a given volume is a fairly different criterion. 

VII. GRANULAR HYDRODYNAMIC THEORY 

A. Derivation 

We take the conserved energy w{s, Sg, p, gi, Uij) of granular media to depend on entropy s, 
granular entropy Sg, density p, momentum density gi, and the elastic strain Uij. Defining the 
conjugate variables as T = dw/ds, Tg = Tg — T = dw/dsg [see Eq (l23l) ]. p — v'^/2 = dw/dp 
[see Eq ([8])], Vi = dw/dgi = gi/ p [see Eq (|5])], 7r,y = —dw/duij, we write 

dw = Tds + Tgdsg + (p - v^/2)dp (61) 

+Vidgi - TTijduij. 
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The equations of motion for the energy and its variables are 

dtw + ViQi = 0, dtp + ViifWi) = 0, (62) 

dtgi + Vj(aij + giVj) = 0, (63) 

dtS + VJi = R/T, dtSg + V,F, = Rg/TG, (64) 

dtUij - Vij - Xij = (65) 

-[{uikVjVk + ViVj/^) + (i ^ j)]. 

The first three equations are conservation laws, with the fluxes Qi and aij as yet unknown, to 
be determined in this section. The next two are the balance equation for the two entropies, 
the form of which are already given, in Eqs (I30II32|33||M1) . Nevertheless, to see that they 
indeed fit the constraints required by energy and momentum conservation, we designate the 
currents as /« = svi — ff, Fi = SgVi — Ff, leaving ff',Ff,R,Rg unspecified. The last is 
the equation of motion for the elastic strain field, as discussed in section |Vl with yi, Xij the 
unknown fluxes to be determined here. Next, we introduce ajj + E^, as 

a,,^(-/ + /.p)5.,-(aJ + Sg) (66) 

+TTij - TCikUjk - TTjkUik, 

where f = Wq —Ts — TgSg, as in Eq (I27|48l) . This is simply a definition of a^ + S^, which 
transfer our task from determining o"jj to finding the new quantity. This simplifies our task, 
notationally, of finding the form of aij, it does not in anyway prejudice it. 

Differentiating the energy, dtw = Tdts + TgdtSg + (/i — v'^ /2)dtp + VidtQi — iTijdtUij, see 
Eq (!6T|) . then inserting Eqs (I62II63II64|[65!) into it, employing relations such as TgdtSg = 
TgRg/Tc + VkSgVkfg - VkifgSgVk), ^6 obtalu 

VjQi = Vi(r/i + TgFi + fxpvi + Vjaij - yj-Kij) (67) 

-R + /f ViT + af^v.j + ViVjTiij + XijTiij + 77;^ 
-Rg + Sgt;,, + F^V.Tg - ^f^ 

This is a useful result, which shows one can rewrite dtW as the divergence of something (first 
line), plus something (second and third line) that vanishes in equilibrium - see section UlI CI 
why VjT, t>y, vTjj, VjTTjj and Tq vanish. We take the first line to yield the energy flux, Qi, 
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and the next two lines to vanish independently, 



Qi = Tfi + TgFi + upvi + Vjffij - yjTiij, (68) 

R = ff'V^T + afjVij + yiVj^l^, + X^.T^i, + -ffl (69) 

Rg = S,%. + i^^V.f, - 7?;^ (70) 

Comparing i?, Rg with Eqs (I32ll34p , the currents are found as 

fF = f^^iT, afj = (vuSij + r]vl + arcij, (71) 

Hi = P^VjTiij, Xij = Pn^j + Pi6ijiTu - avij. 

(It is an assumption to take FfViTg as part of Rg rather than R.) The two terms preceded 
by a contribute iavTjjfjj to i?, respectively, hence cancel each other and are compatible 
with Eq ( |32l) . (More such pairs of terms, mutually canceling or contributing in equal parts, 
are possible. They have been excluded as a simplification. In the language of the Onsager 
force-fiux relation, the above fluxes possess only diagonal elements, with the exception of 
the reactive, off-diagonal terms ~ a.) Defining two relaxation times, 

1^2/?^v^, l = 3Av^(s+^). (72) 

the last of Eqs (I7T1) may be written as 

Xij = A Sij/n - u'^^j/t - avij. (73) 

To ensure permanent elasticity in granular statics, we must in addition require 

Xij -^0 for Tg^O. (74) 

This completes the derivation of GSH. Given ff, Ff,a^,T,l^j,yi,Xij, the structure of all 
currents in the set of equation, Eqs fl62|63|64|[65|) . are known. The question that remains 
is whether these expressions are unique. For simpler hydrodynamic theories, such as for 
isotropic liquid, nematic liquid crystal, or elastic solid, this procedure (frequently referred 
to as the standard procedure) is easily shown to be unique, because one can convince oneself 
that as long as the energy w remains unspecified, there is only one way to write the time 
derivative of the energy dtw as the sum of a divergence and a series of expressions that vanish 
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in equilibrium. In the present case, with two levels of entropy productions, one of which 
controls the switch between permanent and transient elasticity, the hydrodynamic theory is 
singularly intricate, and peripheral ambiguity remains. Nevertheless, displaying energy and 
momentum conservation explicitly, and reducing to liquid and solid hydrodynamics in the 
proper limits, the given set of equations is certainly a viable and consistent theory. 

A more formal way of obtaining the fluxes of Eqs (!7T!) is to define the flux and force 
vectors as Z = {ff, y„ ag, Xij), Zg = {Ff, Eg), Y = (V^T, Vj^ij, Vij, TCij), Yg = {Vifg, Vij). 
And because R = Z -Y, Rg = Zg ■ Yg, the Onsager force-fiux relations are given as 

Z = C-Y, Zg=Cg-Yg. ( 75 ) 

The transport matrices, c and Cg, have positive diagonal elements and off-diagonal ones that 
satisfy the Onsager reciprocity relation. Our example above has only diagonal elements, 
with the single exception of the reactive, off-diagonal terms ~ a. 

B. Results 

Collecting the terms derived above, in section IVII At the equations of GSH, with aij valid 
to lowest order in strain, are 

dtp + Vi{pVi) = 0, (76) 

dfUij = (1 - a)vij - u°j/t - uu 6ij/Ti 
-{u^kVjVk + Vi[/?^Vfc7r,fe/2]) - (i ^ j), (77) 

^ij = (1 ~ '^)^Ji " T^ikUjk - T^jkUik 

+ {fip - f)6,j - (C + Cg)v,A, - (r/ + r]g)v^^, (78) 

Tg[dtSg + Vi{SgVi - KgVifg)] = Rg = 

CA + VgVlvl + ng{V,fgf - ^fl (79) 

T[dts + V^{sv, - «:V.T)] = C4 + VV%vl + ^f^g 

+K{V,Tf + /?^(V,vr,,)2 + /57r° 4 + /^i^. (80) 

Given in terms of the variables: (s, Sg, p, gi, Uij), conjugate variables (T, Tg, p, Vi, ttjj), 
and 11 transport coefficients, (a, r, ti, (, (g, t], rjg, 7, P^, k, Kg), these equations are 
valid irrespective of the functional form of the energy w and the transport coefficients. 

31 



Therefore, they only provide a hydro dynamic structure, a framework into which different 
concrete theories fit. This circumstance, though also true for Newtonian fluids, is not as 
relevant there, because static susceptibilities (such as the compressibility or specific heat) 
and transport coefficients may frequently be approximated as constant. So the structure 
alone already possesses considerable predicting power. This is not true for granular media, 
which typically possess more involved functional dependence - especially concerning the 
Tg —>■ limit, which does not have a counter part in other systems. This is one of the less 
recognized reasons, we believe, underlying the complexity of granular systems. 

In section Wl\ a free energy was proposed that we are confident should be fairly realistic. 
The situation with respect to the 11 transport coefficients are less settled, and in need of 
much future work, though a few limits are clear from the onset: First, a simple, analytic 
way to assure the elastic limit for Tg = and satisfy the requirement of Eq ([71]) is given by 

1/r = AT,, l/n = Aifg, (81) 

which, as we shall see next, gives rise to the same dynamic structure as hypoplasticity. The 
density dependence is more subtle, hence harder and less urgent to determine. However, it 
seems plausible that A, Ai should decrease for growing density, and especially the compres- 
sional relaxation should stop being operative at the random close packing density pcp- To 
account for this, the simplest dependence would be 

Ai~(p-Pcp). (82) 

The coefficient a needs to vanish in the elastic limit, for Tg —>■ 0, and be constant in 
the hypoplastic one, when Tg is moderately large: We have aij = iiij in the elastic regime, 
and aij = (1 — a)TTij + ■ ■ ■ with 1 — a ^ 0.2 in the hypoplastic one, implying sand is much 
softer here - same strain, yet stress is smaller by a factor of about five. The behavior of a 
is probably the result of granular agitation disrupting force chains. They are all intact in 
the elastic limit, making the system comparatively stiff. A finite Tg breaks up the chains, 
and when most of chains are destroyed, the remaining ones become essential in the sense 
that their disruption leads to local collapse, which in turn immediately repair the chains by 
some rearrangement. This is why a saturates and becomes constant. 

Finally, as long as Eq (HO!) holds, the rate independence it entails would prevent the prop- 
agation of sound and elastic waves: Because both the elastic and the plastic part are linear 
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in the velocity, and of the same order in the wave vector g, sound damping is comparable to 
sound velocity, and wave propagation could at most persist for a few periods. We therefore 
expand 7, ?7f, in T^, as 

7 = 7o + 71^9, Vg = ViTg, (83) 

assuming rjg lacks a constant term, because viscous dissipation occurs directly via rj for 
Tg -^ 0, see Eq ( !7T|) . Inserting these expression into Eq (135|) for a quick, qualitative estimate, 
we find Tg ~ VijVij = vl for 70 ^ 71^9, and Tg ~ Vs for 70 <^ liTg- The first regime is 
essentially elastic, because the relaxation term, Uij/r ~ UijTg ~ Mj^f^, is of second order 
and small. This ensures the propagation of sound modes. In the second regime, the same 
term, Uij/t ~ UijTg ~ UijVs-, is of first order and rather more prominent, giving rise to the 
hypoplastic behavior discussed in the next section . 

VIII. THE HYPOPLASTIC REGIME 



Hypoplasticity [20|, a modern, well- verified, yet comparatively simple theory of soil me- 
chanics, models solid dynamics as realistically as the best of the elasto-plastic theories. Its 
starting point is the rate-independent constitutive relation, 

dtaij = Hijke Vke + A^ y t;°4;°- + e vj^, (84) 

where the coefficients Hiji:e,Aij,e are functions of <Jij,p, specified using experimental data 
mainly from triaxial apparatus. (Rate-independence means dtcrij is linearly proportional to 
the magnitude of the velocity, such that the change in stress remains the same for given 
displacement irrespective how fast the change is applied, a well verified observation in both 
the elastic and hypoplastic regime.) Great efforts are invested in finding accurate expressions 
for the coefficients, of which a recent set [20| is e = 1/3, 

Hijke = f {F SikSji + a (Tij(rki/(r^^) , (85) 

Aij = affdF {aij + a") /(t„„, (86) 

where a = 2.76, hs = 1600 MPa, e^ = 0.44ej, Cc = O.SScj, e^^ = exp {au/hg) ' , and 
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GSH, as derived above, reduces to Eq (18^ for a stationary Tg, with Hijik,Aij,e given 
in terms of Mijke = —d'^w/duijduke and four scalars that are combinations of transport 
coefficients. We assume uniformity and stationarity, with especially VjTg, VjTTjj, 9(fj = 0, 
and only include the lowest order terms in the strain Uij. We also take a, rjg, (g as constants, 
and neglect Pt from Eq (135]) . the pressure relevant in granular liquid, assuming Tg is too 
small for the given velocity. It is then quite easy to evaluate dt<7ij employing Eqs (I77II78I) . 

dtaij = (1 - a)dt7Tij = (1 - a)MijkfdtUke = 
(1 - a)MijM[{l - a)vki - uIJt - dkiUmm/ri]. (87) 

Clearly, given Eqs (13511811) . this expression already has the structure of Eq (IM|) that Hy- 
poplasticity postulates. And the coefficients are 

Hijke = (1 - afMijki, e = Cghg, (88) 

A,, = (1 - a)MijkA{T/Ti)^5ki - ul^]\^f^. (89) 

HPM has 43 free parameters (36+6+1 for Hijke, ^ij, e), all functions of the stress and density. 
Expressed as here, the stress and density dependence are essentially determined by Mijke. 



that is a known quantity [6J, l65|. For the four free constants, we take 

(90) 



1- 


-a = 0.22, 


r/n = 0.09, 


Vg 


= 0.33, y^ = 


fm 114 
V7i A 



to be realistic choices, as these numbers yield satisfactory agreement with Hypoplasticity. 
Their significance are: (g/rig = 0.33 implies shear flows are three times as effective in creating 
Tg as compressional flows, t/ti = 0.09 means, plausibly, that the relaxation rate of shear 
stress is ten times higher than that of pressure. The factor (1 — a)^ accounts for an overall 
softening of the static compliance tensor Mijik- Finally, A controls the stress relaxation 
rate for given Tg, and v^tT/Ti how well shear flow excites Tg. Together, X\/%/i =114 
determines the relative weight of plastic versus reactive response: The term in Eq ( 15^ 
preceded by Hijke is the reversible, elastic response, the second term preceded by Ajj comes 
from stress relaxation, is dissipative, irreversible and plastic. For small strain, A,Us -^ 0, 
the elastic part is dominant, \Hijke\ ^ l-^ijl- But \Aij\/\Hijke\ ~ Iti^^l -114/(1 — a) is of order 
unity for \uij\ around 10~^. 
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Although the functions of Eqs fl88ll89p appear rather different from that of Eqs (18511861) . 
the stress-strain increments are quite similar, as the comparison in |66] shows. Moreover, 
the residual discrepancies may be eliminated by discarding the simplifying assumption of 
constant transport coefficients, independent of the stress. This agreement provides valuable 
insights into the physics of Hypoplasticity, showing why it works, what its range of validity 
is, and how it may be generalized. And it conversely also verifies GSH. 

IX. CONCLUSION 

The success of Granular Elasticity, the theory we employ to account for static stress 
distribution in granular media, is mainly due to the fact that the information on the plastic 
strain is quite irrelevant there. This is no longer true in granular dynamics, when the 
system is being deformed - sheared, compressed or tapped. Starting from the working 
hypothesis that granular media are transiently elastic while being deformed, we aim to 
understand the notoriously complex plastic motion by combining two simple and transparent 



elements, elasticity and stress relaxation. In a recently published Letter [66||, we proposed 
a model for granular solids based on this hypothesis. In the present manuscript, we give 
this model a consistent, hydrodynamic framework, compatible with conservation laws and 
thermodynamics. 

The framework is valid for any healthy energy, but is essentially devoid of predictive power 
if the energy is left unspecified. Therefore, an explicit expression for the total, conserved 
energy is given. Encapsulating the key features of static granular media: stress distribution, 
incremental stress-strain relation, minimal and maximal density, the virgin consolidation 
line, the Coulomb yield line and the cap model, this expression should prove realistic enough 
for rendering the specific hydrodynamic theory useful. Much future work is needed to see 
whether further agreement between theory and experiments may be achieved, especially 
concerning cyclic loading, tapping and shear band. 
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APPENDIX A: EQUILIBRIUM CONDITIONS 

First, noting iiijduij = nijdVjUi because iiij is symmetric, we write the energy density 
per unit volume (dropping the subscript of Wq in this section) as 

dw = Tds + /idp — vTjjdVjf/j. (Al) 

Next, we vary the energy J wdV for given entropy J sdV, mass J pdV, and for fixed displace- 
ment at the medium's surface, 6Ui = 0. Taking £1,^2 as constant Lagrange parameters and 
denoting the surface element as dAi, we require the variation of the energy to be extremal. 



6 {w-iis- £2p) dV = 0. (A2) 

Inserting Eq (jAip into (]A2p . we find 

f [TSs + fxSp + TCijSVjUi - ii5s - £25p] dV = 
I [(T - h) 6s + {p- £2) 6p - (V.TT,,) 6U,] dV 

+ (f> TiijSUi dAi = 0, 



where the last term vanishes because 6Ui = at the surface. If 6s, 6p and 6Uj vary 
independently, all three brackets must vanish. And because £±,£2 are constant, T,p also 
need to be. So the conditions for the energy (or entropy) being extremal are 

ViT = 0, V,p = 0, V,7ri, = 0. (A3) 

In granular media for Tg = 0, density and compression are coupled as 

duii = -dp/ p = pdv, (A4) 

and do not vary independently. Simply inserting this relation into the above calculation, 
we find Vi{p + nu/Sp) = to replace the last two conditions of Eq flA3l) . This is not the 
correct result, because we have been varying the energy and its variables above, keeping the 
volume unchanged throughout, with 6Ui = at the surface. But then uu is fixed and cannot 
change with the density p: Consider a one- dimensional medium between x = and x = Xq, 
with Ux{0), Ux{xo) given. Since VjTTjj ~ dl[Ux{xo) — Ux{0)] = 0, the one- dimensional strain 
is Uie = dxUx = {Ux{xq) — Ux{0))/xo and cannot change. 
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To find the proper expression, we may take mass M ratlier than volume V as the constant 
quantity, and vary the density by changing the volume, or the length in the one-dimensional 
case. Holding 5Ui = at the moving surface will then allow Eq (lA4p to hold. Denoting the 
energy, entropy and volume per unit mass, respectively, as e = w / p, a = s/p, v = V/M = 
1/p, and f = w — Ts, the equivalent expression, 

de = Tda - Pxdv - {nij/ p)dVjUi, (A5) 

Pt = -w + Ts + pp = -f + pp, (A6) 

holds. Now we have E = J edM, S = J adM, V = J vdM, where dM = pdV is the 
integrating mass element. Varying the energy for given entropy, volume and requiring it to 
vanish, 6E — iiS — l-^^ = 0, we find 

f [{T - £,) 5a + {Pt - 4) H dM = f (V.-vr,,) SUi dV. 

implying VjT = 0, ViPr = 0, and VjTTjj = 0. These are the same conditions as Eq flA3l) . 
because ViPr = sVjT + pVip. But if Eq (lA4p is implemented, turning Eq (lASp to 

de = Tda - {Pt + 7r«/3)dt; - (7r°./p)dV,t/i, (A7) 

= Tda - p-\PT 5ij + TTij) dVjUi, (A8) 

the equilibrium conditions are altered to become 

ViT = 0, Vi(PT + W3) = 0, V,-7r°. = 0. (A9) 

Clearly, the Cauchy, or total, stress in equilibrium is given as 

aij = PrSij + Tiij, with VjCijj = 0. (AlO) 

Including the gravitational energy p0 in w, with — Vj0 = Gi, the gravitational constant 
pointing downward, we have 

dw = Tds + {p + (f))dp - TTijdVjUi + pd0, (All) 

and find (via the same calculation as above) that /i + is now a constant, implying an 
alteration of the second of Eqs flA3p to 

Vip = Gi, (A12) 
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or equivalently, V-iPt = sVjT + pVi^ = pGi. If Eq f lA4p holds, VjCTjj = is analogously 
changed to 

V,a,, = VjiPrSij + Tiij) = pGi. (A13) 
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